#include "layer.h"

#define SIGN(a) (a>0)?1:0
#define EM_LS  GSL_CONST_MKSA_SPEED_OF_LIGHT // 2.99792458e8	//скорость света
#define M_MAX 80 // number of gn for fi&xi

Layer::Layer()
{
    for(int j=0;j<9; j++) RES[j]=0;
}

Layer::Layer(Structure *ft)
{
    e1=ft->e1;
    e2=ft->e2;
    waveType=ft->waveType;
    Num=ft->difQ[0];
    mI=0; nI=0;
    d=ft->a; a=d/2;   k=ft->k;
    length=2*M_PI/k;
    R=ft->R; teta=ft->teta; fi=ft->fi; Z=ft->Z; Dist=ft->Dist;
    for (int j=0; j<17; j++) { Tm[j]=0; Her[j]=0; }
    for(int j=0;j<9; j++) RES[j]=0;
    roots();
    coeffs();
}

gsl_complex Layer::pole(double x, double z)
{
    gsl_complex EH, gslBuf;
    double Rs,rho,Phi;
    //R=z*(1+pow(M_PI/length*pow(a,2)/z,2));
    Rs=z*M_PI*length/(pow(length*z,2)+pow(M_PI*pow(a,2),2));
    //rho=a*(1+pow(length/M_PI/pow(a,1)*z,2));    //before demistake
    rho=a*sqrt(1+pow(length/M_PI/pow(a,2)*z,2));
    Phi=atan(length/M_PI/pow(a,2)*z);
    //getHermitOdd(sqrt(2)*x/rho);
    getHermitEven(sqrt(2)*x/rho);
    GSL_SET_COMPLEX(&EH,0,0);
    for (int j=0; j<Num; j++)
    {
        //GSL_SET_COMPLEX(&gslBuf,-pow(x/rho,2),k*z-Phi*(j+1)+pow(x,2)*Rs);    //before demistake
        GSL_SET_COMPLEX(&gslBuf,-pow(x/rho,2),k*z-(j+1)*Phi+pow(x,2)*Rs);
        gslBuf=gsl_complex_mul_real(gsl_complex_exp(gslBuf),Tm[j]*a/rho*Her[j]);
        EH=gsl_complex_add(EH,gslBuf);
    }
    return EH;
}

void Layer::roots()
{
    int j=0;
    double alf=0,xa,xb,det,V;
    __complex__ double d1,d2;
    V=pow(k*d,2)*(e1-e2);
    det=0.1;
    while (j<1)
    {
        xa=alf;	xb=xa+det;
        d1=disper(xa);
        d1-=V;
        d2=disper(xb);
        d2-=V;
        if (GSL_SIGN(__real__(d1))!=GSL_SIGN(__real__(d2)))
        {
            for (int q=0;q<16;q++)
            {
                alf=(xa+xb)/2;
                d1=disper(xa)-V;
                d2=disper(alf)-V;
                if (GSL_SIGN(__real__(d1))!=GSL_SIGN(__real__(d2)))	xb=alf;
                else	xa=alf;
            }
            u[j]=alf;
            switch(waveType)
            {
            case 1 : // E-type
                v[j]=u[j]*tan(u[j]/2)*e2/e1;
                break;
            case 2 : // H-type
                v[j]=u[j]*tan(u[j]/2);
                break;
            }
            j++;
        }//j++;
        alf=alf+det;
    }
}

double Layer::disper(double x)
{
    double val;
    switch(waveType)
    {
    case 1 : // E-type
        val=pow(x,2)+pow(x*tan(x/2)*e2/e1,2);
        break;
    case 2 : // H-type
        val=tan(x/2);
        val*=val;
        val+=1;
        val*=pow(x,2);
        break;
    }
    return val;
}

void Layer::coeffs()
{
    int N=100;
    double Buf[17],x,dx,e,buf;
    for (int j=0; j<=Num; j++) Buf[j]=0;
    dx=a/N*2;
    x=a;
    for (int i=0; i<N/2; i++)
    {
        getHermitEven(sqrt(2)*x/a);
        e=exp(-pow(x/a,2));
        for (int j=0; j<=Num; j++)
        {
            if (remainder(j,2)==0)
            {
                buf=cos(u[0]*x/(2*a));
                Buf[j]+=buf*Her[j]*e;
            }
        }
        x-=dx;
    }
    for (int j=0; j<=Num; j++) Tm[j]=Buf[j]*dx;
    for (int j=0; j<=Num; j++) Buf[j]=0;
    dx=a/N*5;
    x=a;
    for (int j=0; j<N; j++)
    {
        getHermitEven(sqrt(2)*x/a);
        for (int j=0; j<=Num; j++)
        {
            if (remainder(j,2)==0)
            {
                e=exp(-v[0]/(2*a)*(x-a)-pow(x/a,2));
                Buf[j]+=Her[j]*e;
            }
        }
        x+=dx;
    }
    for (int j=0; j<=Num; j++)
    {
        Tm[j]+=cos(u[0]/2)*Buf[j]*dx;
        Tm[j]/=sqrt(M_PI)*pow(2,j)*gsl_sf_fact(j)/2;
    }
}

//void Layer::getHermitOdd(double x)
//{
//    double z[10];
//    z[0]=1;
//    for (int i=1; i<10; i++)
//    {
//        if (i<Num) z[i]=z[i-1]*x;
//        else z[i]=0;
//    }
//    Her[1]=2*z[1];
//    Her[3]=8*z[3]-12*z[1];
//    Her[5]=32*z[5]-160*z[3]+120*z[1];
//    Her[7]=128*z[7]-1344*z[5]+3360*z[3]-1680*z[1];
//}

void Layer::getHermitEven(double x)
{
    double z[17];
    z[0]=1;
    for (int i=1; i<=16; i++)
    {
        if (i<=Num) z[i]=z[i-1]*x;
        else z[i]=0;
    }
    Her[0]=1;
    Her[2]=4*z[2]-2;
    Her[4]=16*z[4]-48*z[2]+12;
    Her[6]=64*z[6]-480*z[4]+720*z[2]-120;
    Her[8]=256*z[8]-3584*z[6]+13440*z[4]-13440*z[2]+1680;
    Her[10]=1024*z[10]-23040*z[8]+161280*z[6]-403200*z[4]+302400*z[2]-30240;
    Her[12]=4096*z[12]-135168*z[10]+1520640*z[8]-7096320*z[6]+13305600*z[4]-7983360*z[2]+665280;
    Her[14]=16384*z[14]-745472*z[12]+12300288*z[10]-92252160*z[8]+322882560*z[6]-484323840*z[4]+242161920*z[2]-17297280;
    Her[16]=65536*z[16]-3932160*z[14]+89456640*z[12]-984023040*z[10]+5.535129600*1e9*z[8]-1.5498362880*1e10*z[6]+1.9372953600*1e10*z[4]-8.302694400*1e9*z[2]+518918400;
}

__complex__ double Layer::compSqrt(__complex__ double val)
{
    gsl_complex gslBuf;
    GSL_SET_COMPLEX(&gslBuf,__real__(val),__imag__(val));
    gslBuf=gsl_complex_sqrt(gslBuf);
    val=GSL_REAL(gslBuf)+1I*GSL_IMAG(gslBuf);
    return val;
}
